home *** CD-ROM | disk | FTP | other *** search
/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / cpp_libs / lin_alg.lha / lin_alg / vvector.cc < prev    next >
C/C++ Source or Header  |  1993-08-08  |  11KB  |  378 lines
  1. // This may look like C code, but it is really -*- C++ -*-
  2. /*
  3.  ************************************************************************
  4.  *
  5.  *            Verify Operations on Vectors
  6.  *
  7.  ************************************************************************
  8.  */
  9.  
  10. #include "LinAlg.h"
  11. #include <math.h>
  12. #include <builtin.h>
  13. #include <ostream.h>
  14.  
  15. /*
  16.  *------------------------------------------------------------------------
  17.  *            Service validation functions
  18.  */
  19.  
  20. static void verify_element_value(const Vector& v,const REAL val)
  21. {
  22.   register imax = 0;
  23.   register double max_dev = 0;
  24.   register int i;
  25.   for(i=v.q_lwb(); i<=v.q_upb(); i++)
  26.   {
  27.     register double dev = abs(v(i)-val);
  28.     if( dev >= max_dev )
  29.       imax = i, max_dev = dev;
  30.   }
  31.  
  32.   if( max_dev == 0 )
  33.     return;
  34.   else if( max_dev < 1e-5 )
  35.     message("Element #%d with value %g differs the most from what\n"
  36.         "was expected, %g, though the deviation %g is small\n",
  37.         imax,v(imax),val,max_dev);
  38.   else
  39.     _error("A significant difference from the expected value %g\n"
  40.        "encountered for element #%d with value %g",
  41.        val,imax,v(imax));
  42. }
  43.  
  44. static void verify_vector_identity(const Vector& v1, const Vector& v2)
  45. {
  46.   register imax = 0;
  47.   register double max_dev = 0;
  48.   register int i;
  49.   are_compatible(v1,v2);
  50.   for(i=v1.q_lwb(); i<=v1.q_upb(); i++)
  51.   {
  52.     register double dev = abs(v1(i)-v2(i));
  53.     if( dev >= max_dev )
  54.       imax = i, max_dev = dev;
  55.   }
  56.  
  57.   if( max_dev == 0 )
  58.     return;
  59.   if( max_dev < 1e-5 )
  60.     message("Two #%d elements of the vectors with values %g and %g\n"
  61.         "differ the most, though the deviation %g is small\n",
  62.         imax,v1(imax),v2(imax),max_dev);
  63.   else
  64.     _error("A significant difference between the vectors encountered\n"
  65.        "at #%d element, with values %g and %g",
  66.        imax,v1(imax),v2(imax));
  67. }
  68.  
  69. /*
  70.  *------------------------------------------------------------------------
  71.  *       Test allocation functions and compatibility check
  72.  */
  73.  
  74. static void test_allocation(void)
  75. {
  76.  
  77.   cout << "\n\n---> Test allocation and compatibility check\n";
  78.  
  79.   Vector v1(20);
  80.   Vector v2(1,20);
  81.   Vector v3(0,19);
  82.   Vector v4(v1);
  83.   cout << "The following vector have been allocated\n";
  84.   v1.info(), v2.info(), v3.info(), v4.info();
  85.  
  86.   cout << "\nStatus information reported for vector v3:\n";
  87.   cout << "  Lower bound ... " << v3.q_lwb() << "\n";
  88.   cout << "  Upper bound ... " << v3.q_upb() << "\n";
  89.   cout << "  No. of elements " << v3.q_no_elems() << "\n";
  90.   cout << "  Name " << v3.q_name() << "\n";
  91.  
  92.   cout << "\nCheck vectors 1 & 2 for compatibility\n";
  93.   are_compatible(v1,v2);
  94.  
  95.   cout << "Check vectors 1 & 4 for compatibility\n";
  96.   are_compatible(v1,v4);
  97.  
  98. #if 0
  99.   cout << "Check vectors 1 & 3 for compatibility\n";
  100.   are_compatible(v1,v3);
  101. #endif
  102.  
  103.   cout << "v2 has to be compatible with v3 after resizing to v3\n";
  104.   v2.resize_to(v3);
  105.   are_compatible(v2,v3);
  106.  
  107.   Vector v5(v1.q_upb()+5); v5.set_name("Vector v5");
  108.   cout << "v1 has to be compatible with v5 after resizing to v5.upb\n";
  109.   v5.info();
  110.   v1.resize_to(v5.q_no_elems());
  111.   are_compatible(v1,v5);
  112.  
  113.   cout << "\nDone\n";
  114. }
  115.  
  116. /*
  117.  *------------------------------------------------------------------------
  118.  *             Test uniform element operations
  119.  */
  120.  
  121. static void test_element_op(const int vsize)
  122. {
  123.   const double pattern = PI;
  124.   register int i;
  125.  
  126.   cout << "\n\n---> Test operations that treat each element uniformly\n";
  127.  
  128.   Vector v(-1,vsize-2);
  129.   Vector v1(v);
  130.  
  131.   cout << "Writing zeros to v...\n";
  132.   for(i=v.q_lwb(); i<=v.q_upb(); i++)
  133.     v(i) = 0;
  134.   verify_element_value(v,0);
  135.  
  136.   cout << "Clearing v1 ...\n";
  137.   v1.clear();
  138.   verify_element_value(v1,0);
  139.  
  140.   cout << "Comparing v1 with 0 ...\n";
  141.   assure(v1 == 0, "v1 is not zero!");
  142.  
  143.   cout << "Writing a pattern " << pattern << " by assigning to v(i)...\n";
  144.   for(i=v.q_lwb(); i<=v.q_upb(); i++)
  145.     v(i) = pattern;
  146.   verify_element_value(v,pattern);
  147.  
  148.   cout << "Writing the pattern by assigning to v1 as a whole ...\n";
  149.   v1 = pattern;
  150.   verify_element_value(v1,pattern);
  151.  
  152.   cout << "Comparing v and v1 ...\n";
  153.   assure(v == v1, "v and v1 are unexpectedly different!");
  154.   cout << "Comparing (v=0) and v1 ...\n";
  155.   assert(!(v.clear() == v1));
  156.  
  157.   cout << "\nClear v and add the pattern\n";
  158.   v.clear();
  159.   v += pattern;
  160.   verify_element_value(v,pattern);
  161.   cout << "   add the doubled pattern with the negative sign\n";
  162.   v += -2*pattern;
  163.   verify_element_value(v,-pattern);
  164.   cout << "   subtract the trippled pattern with the negative sign\n";
  165.   v -= -3*pattern;
  166.   verify_element_value(v,2*pattern);
  167.  
  168.   cout << "\nVerify comparison operations\n";
  169.   v = pattern;
  170.   cout << "   (v=pattern)>0\n";
  171.   assure( v > 0, "Unexpectedly, v is not greater than 0");
  172.   cout << "   (v=pattern)>-pattern\n";
  173.   assert( v > -pattern );
  174.   cout << "   (v=-pattern)<-pattern/2\n";
  175.   v -= 2*pattern;
  176.   assert( v  < -pattern/2 );
  177.  
  178.   cout << "\nAssign 2*pattern to v by repeating additions\n";
  179.   v = 0; v += pattern; v += pattern;
  180.   cout << "Assign 2*pattern to v1 by multiplying by two \n";
  181.   v1 = pattern; v1 *= 2;
  182.   verify_element_value(v1,2*pattern);
  183.   assert( v == v1 );
  184.   cout << "Multiply v1 by one half returning it to the 1*pattern\n";
  185.   v1 *= 1/2.;
  186.   verify_element_value(v1,pattern);
  187.  
  188.   cout << "\nAssign -pattern to v and v1\n";
  189.   v.clear(); v -= pattern; v1 = -pattern;
  190.   verify_element_value(v,-pattern);
  191.   assert( v == v1 );
  192.   cout << "v = sqrt(sqr(v)); v1 = abs(v1); Now v and v1 have to be the same\n";
  193.   v.sqr();
  194.   verify_element_value(v,pattern*pattern);
  195.   v.sqrt();
  196.   verify_element_value(v,pattern);
  197.   v1.abs();
  198.   verify_element_value(v1,pattern);
  199.   assert( v == v1 );
  200.  
  201.   cout << "\nCheck out to see that sin^2(x) + cos^2(x) = 1\n";
  202.   for(i=v.q_lwb(); i<=v.q_upb(); i++)
  203.     v(i) = 2*PI/v.q_no_elems() * i;
  204.   v1 = v;
  205.   v.user_function(sin);
  206.   v1.user_function(cos);
  207.   v.sqr();
  208.   v1.sqr();
  209.   v += v1;
  210.   verify_element_value(v,1);
  211.  
  212.   cout << "\nVerify constructor with initialization\n";
  213.   Vector vi(1,5,0.0,1.0,2.0,3.0,4.0,"END");
  214.   Vector vit(5);
  215.   for(i=vit.q_lwb(); i<=vit.q_upb(); i++)
  216.     vit(i) = i-1;
  217.   verify_vector_identity(vi,vit);
  218.   
  219.   cout << "\nDone\n\n";
  220. }
  221.  
  222. /*
  223.  *------------------------------------------------------------------------
  224.  *            Test binary vector operations
  225.  */
  226.  
  227. static void test_binary_op(const int vsize)
  228. {
  229.   const double pattern = PI;
  230.   register int i;
  231.  
  232.   cout << "\n---> Test Binary Vector operations\n";
  233.  
  234.   Vector v(2,vsize+1);
  235.   Vector v1(v);
  236.   Vector vp(v);
  237.  
  238.   for(i=v.q_lwb(); i<=v.q_upb(); i++)
  239.     vp(i) = (i-v.q_no_elems()/2.)*pattern;
  240.   
  241.  
  242.   cout << "\nVerify assignment of a vector to the vector\n";
  243.   v = pattern;
  244.   v1.clear();
  245.   v1 = v;
  246.   verify_element_value(v1,pattern);
  247.   assert( v1 == v );
  248.  
  249.   cout << "\nAdding one vector to itself, uniform pattern " << pattern << "\n";
  250.   v.clear(); v = pattern;
  251.   v1 = v; v1 += v1;
  252.   verify_element_value(v1,2*pattern);
  253.   cout << "  subtracting two vectors ...\n";
  254.   v1 -= v;
  255.   verify_element_value(v1,pattern);
  256.   cout << "  subtracting the vector from itself\n";
  257.   v1 -= v1;
  258.   verify_element_value(v1,0);
  259.   cout << "  adding two vectors together\n";
  260.   v1 += v;
  261.   verify_element_value(v1,pattern);
  262.  
  263.   cout << "\nArithmetic operations on vectors with not the same elements\n";
  264.   cout << "   adding vp to the zero vector...\n";
  265.   v.clear(); v += vp;
  266. //  verify_vector_identity(v,vp);
  267.   assert( v == vp );
  268.   v1 = v;
  269.   cout << "   making v = 3*vp and v1 = 3*vp, via add() and succesive mult\n";
  270.   add(v,2,vp);
  271.   v1 += v1; v1 += vp;
  272.   verify_vector_identity(v,v1);
  273.   cout << "   clear both v and v1, by subtracting from itself and via add()\n";
  274.   v1 -= v1;
  275.   add(v,-3,vp);
  276.   verify_vector_identity(v,v1);
  277.  
  278.   cout << "\nTesting element-by-element multiplications and divisions\n";
  279.   cout << "   squaring each element with sqr() and via multiplication\n";
  280.   v = vp; v1 = vp;
  281.   v.sqr();
  282.   element_mult(v1,v1);
  283.   verify_vector_identity(v,v1);
  284.   cout << "   compare (v = pattern^2)/pattern with pattern\n";
  285.   v = pattern; v1 = pattern;
  286.   v.sqr();
  287.   element_div(v,v1);
  288.   verify_element_value(v,pattern);
  289.   compare(v1,v,"Original vector and vector after squaring and dividing");
  290.  
  291.   cout << "\nDone\n";
  292. }
  293.  
  294. /*
  295.  *------------------------------------------------------------------------
  296.  *            Verify the norm calculation
  297.  */
  298.  
  299. static void test_norms(const int vsize)
  300. {
  301.   cout << "\n---> Verify norm calculations\n";
  302.  
  303.   const double pattern = 10.25;
  304.  
  305.   if( vsize % 2 == 1 )
  306.     _error("Sorry, size of the vector to test must be even for this test\n");
  307.  
  308.   Vector v(vsize);
  309.   Vector v1(v);
  310.  
  311.   cout << "\nAssign " << pattern << " to all the elements and check norms\n";
  312.   v = pattern;
  313.   cout << "  1. norm should be pattern*no_elems\n";
  314.   assert( v.norm_1() == pattern*v.q_no_elems() );
  315.   cout << "  Square of the 2. norm has got to be pattern^2 * no_elems\n";
  316.   assert( v.norm_2_sqr() == sqr(pattern)*v.q_no_elems() );
  317.   cout << "  Inf norm should be pattern itself\n";
  318.   assert( v.norm_inf() == pattern );
  319.   cout << "  Scalar product of vector by itself is the sqr(2. vector norm)\n";
  320.   assert( v.norm_2_sqr() == v * v );
  321.  
  322.   double ap_step = 1;
  323.   double ap_a0   = -pattern;
  324.   int n = v.q_no_elems();
  325.   cout << "\nAssign the arithm progression with 1. term " << ap_a0 <<
  326.           "\nand the difference " << ap_step << "\n";
  327.   register int i;
  328.   for(i=v.q_lwb(); i<=v.q_upb(); i++)
  329.     v(i) = (i-v.q_lwb())*ap_step + ap_a0;
  330.   int l = min(max((int)ceil(-ap_a0/ap_step),0),n);
  331.   double norm = (2*ap_a0 + (l+n-1)*ap_step)/2*(n-l) +
  332.     (-2*ap_a0-(l-1)*ap_step)/2*l;
  333.   cout << "  1. norm should be " << norm << "\n";
  334.   assert( v.norm_1() == norm );
  335.   norm = n*( sqr(ap_a0) + ap_a0*ap_step*(n-1) + sqr(ap_step)*(n-1)*(2*n-1)/6);
  336.   cout << "  Square of the 2. norm has got to be "
  337.           "n*[ a0^2 + a0*q*(n-1) + q^2/6*(n-1)*(2n-1) ], or " << norm << "\n";
  338.   assert( v.norm_2_sqr() == norm );
  339.   norm = max(abs(v(v.q_lwb())),abs(v(v.q_upb())));
  340.   cout << "  Inf norm should be max(abs(a0),abs(a0+(n-1)*q)), ie " << norm <<
  341.           "\n";
  342.   assert( v.norm_inf() == norm );
  343.   cout << "  Scalar product of vector by itself is the sqr(2. vector norm)\n";
  344.   assert( v.norm_2_sqr() == v * v );
  345.  
  346.   v1.clear();
  347.   compare(v,v1,"Compare the vector v with a zero vector");
  348.  
  349.   cout << "\nConstruct v1 to be orthogonal to v as v(n), -v(n-1), v(n-2)...\n";
  350.   for(i=0; i<v1.q_no_elems(); i++)
  351.     v1(i+v1.q_lwb()) = v(v.q_upb()-i) * ( i % 2 == 1 ? -1 : 1 );
  352.   cout << "||v1|| has got to be equal ||v|| regardless of the norm def\n";
  353.   assert( v1.norm_1() == v.norm_1() );
  354.   assert( v1.norm_2_sqr() == v.norm_2_sqr() );
  355.   assert( v1.norm_inf() == v.norm_inf() );
  356.   cout << "But the scalar product has to be zero\n";
  357.   assert( v1 * v == 0 );
  358.  
  359.   cout << "\nDone\n";
  360. }
  361.  
  362. /*
  363.  *------------------------------------------------------------------------
  364.  *                Root module
  365.  */
  366.  
  367. main()
  368. {
  369.   cout << "\n\n" << _Minuses << 
  370.           "\n\t\tVerify Operations on Vectors\n";
  371.  
  372.   test_allocation();
  373.   test_element_op(20);
  374.   test_binary_op(20);
  375.   test_norms(20);
  376. }
  377.  
  378.